### Alizade, Dancygier, Ditlmann 
### "National Penalties Reversed"
### Replication Code 
### Table 6
### For questions, contact jalizade@princeton.edu

# empty environment
rm(list = ls())

#setwd("")
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")

# load necessary packages
library(readstata13)
library(boot)

# load data set
dat <- read.dta13("data_experimental.dta")

# convert treatment variable from second experiment to factor
dat$e2_treat[dat$e2_treat=="NA"] <- NA
dat$e2_treat <- factor(dat$e2_treat, levels=c("Vote Intention (Canadian)", "Control (Turkish)", "Vote Intention (Turkish)", "Integration Problems (Turkish)"))


### table ###

## N
N <- aggregate(dat$e2_response ~ dat$turk_pop_dummy + dat$e2_treat, FUN=function(x) length(x))
names(N) <- c("turk_pop", "treat_cond", "n")
N$turk_pop <- rep(c("Small", "Large"), 4)

## response rate
resp_rate <- aggregate(dat$e2_response ~ dat$turk_pop_dummy + dat$e2_treat, FUN=function(x) mean(x))[,3]
resp_rate <- round(resp_rate*100, 2)

## ratio
ratio <- c(mean(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)" & dat$turk_pop_dummy==1], na.rm=T) / mean(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)" & dat$turk_pop_dummy==0], na.rm=T), NA,
           mean(dat$e2_response[dat$e2_treat=="Control (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) / mean(dat$e2_response[dat$e2_treat=="Control (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA,
           mean(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) / mean(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA,
           mean(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) / mean(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA)
ratio <- round(ratio, 2)

## bootstrap CI for ratio

set.seed(9447)

# canadian vote intention
ratio_canvote <- function(dat, i) {
  dat2 <- dat[i,]
  (mean(dat2$e2_response[dat2$turk_pop_dummy==1 & dat2$e2_treat=="Vote Intention (Canadian)"], na.rm=T) / (mean(dat2$e2_response[dat2$turk_pop_dummy==0 & dat2$e2_treat=="Vote Intention (Canadian)"], na.rm=T)))
}
boot_ratio_canvote <- boot(data=dat, statistic=ratio_canvote, R=10000)
ratio_ci_canvote <- round(boot.ci(boot_ratio_canvote, conf = 0.95, type="perc")$percent[4:5], 2)

# turkish control
ratio_turkcontrol <- function(dat, i) {
  dat2 <- dat[i,]
  (mean(dat2$e2_response[dat2$turk_pop_dummy==1 & dat2$e2_treat=="Control (Turkish)"], na.rm=T) / (mean(dat2$e2_response[dat2$turk_pop_dummy==0 & dat2$e2_treat=="Control (Turkish)"], na.rm=T)))
}
boot_ratio_turkcontrol <- boot(data=dat, statistic=ratio_turkcontrol, R=10000)
ratio_ci_turkcontrol <- round(boot.ci(boot_ratio_turkcontrol, conf = 0.95, type="perc")$percent[4:5], 2)

# turkish vote intention
ratio_turkvote <- function(dat, i) {
  dat2 <- dat[i,]
  (mean(dat2$e2_response[dat2$turk_pop_dummy==1 & dat2$e2_treat=="Vote Intention (Turkish)"], na.rm=T) / (mean(dat2$e2_response[dat2$turk_pop_dummy==0 & dat2$e2_treat=="Vote Intention (Turkish)"], na.rm=T)))
}
boot_ratio_turkvote <- boot(data=dat, statistic=ratio_turkvote, R=10000)
ratio_ci_turkvote <- round(boot.ci(boot_ratio_turkvote, conf = 0.95, type="perc")$percent[4:5], 2)

# turkish integration problems
ratio_turkint <- function(dat, i) {
  dat2 <- dat[i,]
  (mean(dat2$e2_response[dat2$turk_pop_dummy==1 & dat2$e2_treat=="Integration Problems (Turkish)"], na.rm=T) / (mean(dat2$e2_response[dat2$turk_pop_dummy==0 & dat2$e2_treat=="Integration Problems (Turkish)"], na.rm=T)))
}
boot_ratio_turkint <- boot(data=dat, statistic=ratio_turkint, R=10000)
ratio_ci_turkint <- round(boot.ci(boot_ratio_turkint, conf = 0.95, type="perc")$percent[4:5], 2)

# combine CIs into one data frame
ratio_ci <- rbind.data.frame(ratio_ci_canvote, c(NA, NA), ratio_ci_turkcontrol, c(NA, NA), ratio_ci_turkvote, c(NA, NA), ratio_ci_turkint, c(NA, NA))
names(ratio_ci) <- c("ratio_CI_lower", "ratio_CI_upper")

## difference in means
diff <- c(mean(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)" & dat$turk_pop_dummy==1], na.rm=T) - mean(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)" & dat$turk_pop_dummy==0], na.rm=T), NA,
          mean(dat$e2_response[dat$e2_treat=="Control (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) - mean(dat$e2_response[dat$e2_treat=="Control (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA,
          mean(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) - mean(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA,
          mean(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)" & dat$turk_pop_dummy==1], na.rm=T) - mean(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)" & dat$turk_pop_dummy==0], na.rm=T), NA)
diff <- round(diff*100, 2)


## standard error of difference
# perform t-test first
diff_t <- c(t.test(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)"]~dat$turk_pop_dummy[dat$e2_treat=="Vote Intention (Canadian)"])$statistic, NA,
            t.test(dat$e2_response[dat$e2_treat=="Control (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Control (Turkish)"])$statistic, NA,
            t.test(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Vote Intention (Turkish)"])$statistic, NA,
            t.test(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Integration Problems (Turkish)"])$statistic, NA)

# to obtain the standard error, divide the mean difference by the t-statistic
diff_se <- round(abs(diff / diff_t), 3)

## p-value
diff_p <- c(t.test(dat$e2_response[dat$e2_treat=="Vote Intention (Canadian)"]~dat$turk_pop_dummy[dat$e2_treat=="Vote Intention (Canadian)"])$p.value, NA,
            t.test(dat$e2_response[dat$e2_treat=="Control (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Control (Turkish)"])$p.value, NA,
            t.test(dat$e2_response[dat$e2_treat=="Vote Intention (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Vote Intention (Turkish)"])$p.value, NA,
            t.test(dat$e2_response[dat$e2_treat=="Integration Problems (Turkish)"]~dat$turk_pop_dummy[dat$e2_treat=="Integration Problems (Turkish)"])$p.value, NA)
diff_p <- round(diff_p, 3)

## create table
table <- data.frame(N, resp_rate, ratio, ratio_ci, diff, diff_se, diff_p)
table